data Sim;
   keep Specialty Diagnosis Physician Payout;
   array PhysicianEffect{3000};
   array PhysicianDiagnosis{3000};
   array PhysicianSpecialty{3000};
   do i = 1 to dim(PhysicianEffect);
      PhysicianEffect{i}  = sqrt(4.0)*rannor(1234);
      PhysicianDiagnosis{i}    = 1 + int(100*ranuni(1234));
      PhysicianSpecialty{i} = 1 + int(5*ranuni(1234));
   end;
   do i = 1 to 40000;
      Physician = 1 + int(3000*ranuni(1234));
      Specialty = PhysicianSpecialty{Physician};
      Diagnosis = PhysicianDiagnosis{Physician};
      Payout    = (100 + Specialty + Diagnosis/9.9 + PhysicianEffect{Physician}
                  + sqrt(7.8)*rannor(1235))*3;

      format payout dollar9.2;
      output;
   end;
run;
ods graphics off;

proc hpmixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty;
   random Physician/cl;
   test Specialty Diagnosis*Specialty;
run;

proc mixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty;
   random Physician/cl;
run;
/* THIS IS SLOW */
proc mixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty;
   random int/subject=Physician cl;
run;

/* USING EFFICIENT CODING BEST PRACTICES */
proc sort data = sim out = sim2;
   by physician;
run;
proc mixed data=Sim2;
   class Specialty Diagnosis;
   model Payout = Specialty Diagnosis*Specialty/ddfm=residual;
   random int/subject=Physician cl;
run;

proc hpmixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty/S;
   random Physician/cl;
   test Specialty Diagnosis*Specialty;
run;
proc mixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty/S;
   random Physician/cl;
run;



proc hpmixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty/S;
   random Physician/cl;
   test Specialty Diagnosis*Specialty;
   ods output SolutionR=eblup1;
   ods output covparms=a;
run;
proc mixed data=Sim;
   class Specialty Diagnosis Physician;
   model Payout = Specialty Diagnosis*Specialty/S;
   random Physician/cl;
   parms /pdata=a hold=1,2 noiter;
   lsmeans specialty / pdiff=all adjust=tukey;
run;

ods graphics;